home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 4 / Amiga Tools 4.iso / grafix / raytracing / raylab / source / algebra.c next >
C/C++ Source or Header  |  1996-01-04  |  4KB  |  236 lines
  1. /*
  2.     name:    algebra.c
  3.  
  4.     Linear algebra
  5.     --------------
  6.  
  7.     These functions simplify vector algebra etc.
  8.  
  9. */
  10.  
  11. #include  "defs.h"
  12.  
  13.  
  14.  
  15. void CreateVector(VECTOR *v, double vx, double vy, double vz)
  16. {
  17.     v->x=vx;
  18.     v->y=vy;
  19.     v->z=vz;
  20. }
  21.  
  22.  
  23.  
  24. /* v2 = v1 */
  25.  
  26. void CopyVector(VECTOR *v2, VECTOR *v1)
  27. {
  28.     v2->x=v1->x;
  29.     v2->y=v1->y;
  30.     v2->z=v1->z;
  31. }
  32.  
  33.  
  34.  
  35. /* v3 = v1 + v2 */
  36.  
  37. void AddVector(VECTOR *v3, VECTOR *v1, VECTOR *v2)
  38. {
  39.     v3->x=(v1->x)+(v2->x);
  40.     v3->y=(v1->y)+(v2->y);
  41.     v3->z=(v1->z)+(v2->z);
  42. }
  43.  
  44.  
  45.  
  46. /* v3 = v1 - v2 */
  47.  
  48. void SubVector(VECTOR *v3, VECTOR *v1, VECTOR *v2)
  49. {
  50.     v3->x=(v1->x)-(v2->x);
  51.     v3->y=(v1->y)-(v2->y);
  52.     v3->z=(v1->z)-(v2->z);
  53. }
  54.  
  55.  
  56.  
  57. /* v2 = -v1 */
  58.  
  59. void NegVector(VECTOR *v2, VECTOR *v1)
  60. {
  61.     v2->x=-(v1->x);
  62.     v2->y=-(v1->y);
  63.     v2->z=-(v1->z);
  64. }
  65.  
  66.  
  67.  
  68. /* v3 = v1 x v2 */
  69.  
  70. void CrossProduct(VECTOR *v3, VECTOR *v1, VECTOR *v2)
  71. {
  72.     v3->x=(v1->y)*(v2->z)-(v1->z)*(v2->y);
  73.     v3->y=(v1->z)*(v2->x)-(v1->x)*(v2->z);
  74.     v3->z=(v1->x)*(v2->y)-(v1->y)*(v2->x);
  75. }
  76.  
  77.  
  78.  
  79. /* Returns  v1 · v2 */
  80.  
  81. double DotProduct(VECTOR *v1, VECTOR *v2)
  82. {
  83.     return((v1->x)*(v2->x)+(v1->y)*(v2->y)+(v1->z)*(v2->z));
  84. }
  85.  
  86.  
  87.  
  88. /* v2 = t * v1 */
  89.  
  90. void ScaleVector(VECTOR *v2, double t, VECTOR *v1)
  91. {
  92.     v2->x=t*(v1->x);
  93.     v2->y=t*(v1->y);
  94.     v2->z=t*(v1->z);
  95. }
  96.  
  97.  
  98.  
  99. /* Returns |v| */
  100.  
  101. double VectorLength(VECTOR *v)
  102. {
  103.     return(sqrt((v->x)*(v->x)+(v->y)*(v->y)+(v->z)*(v->z)));
  104. }
  105.  
  106.  
  107.  
  108. /* Returns the (smallest) angle between v1 and v2 */
  109.  
  110. double VectorsAngle(VECTOR *v1, VECTOR *v2)
  111. {
  112.     return(acos(DotProduct(v1,v2)/(VectorLength(v1)*VectorLength(v2))));
  113. }
  114.  
  115.  
  116.  
  117. void CreatePoint(POINT *p, double px, double py, double pz)
  118. {
  119.     p->x=px;
  120.     p->y=py;
  121.     p->z=pz;
  122. }
  123.  
  124.  
  125.  
  126. /* p2 = p1 */
  127.  
  128. void CopyPoint(POINT *p2, POINT *p1)
  129. {
  130.     p2->x=p1->x;
  131.     p2->y=p1->y;
  132.     p2->z=p1->z;
  133. }
  134.  
  135.  
  136.  
  137. /* This routine creates a line from two given points */
  138.  
  139. void MakeLine(LINE *l, POINT *p1, POINT *p2)
  140. {
  141.     CopyPoint(&(l->Origin),p1);
  142.     CreateVector(&(l->Direction),(p2->x)-(p1->x),(p2->y)-(p1->y),(p2->z)-(p1->z));
  143. }
  144.  
  145.  
  146.  
  147. /* This routine creates a reflected vector v2, given a vector v1 and a normal n */
  148.  
  149. void ReflectVector(VECTOR *v2, VECTOR *v1, VECTOR *n)
  150. {
  151.     double    k,a;
  152.     VECTOR    vtemp,vtemp2;
  153.  
  154.     CopyVector(&vtemp2,n);
  155.     a=PI-VectorsAngle(v1,n);
  156.     if((a<0.0)&&(a>-EPSILON)) a+=EPSILON;
  157.     if(a>PID2) {
  158.         NegVector(&vtemp2,n);
  159.         a=PI-a;
  160.     }
  161.     k=VectorLength(&vtemp2)/(2*cos(a)*VectorLength(v1));
  162.     ScaleVector(&vtemp,k,v1);
  163.     AddVector(v2,&vtemp,&vtemp2);
  164. /*
  165.     if(k<0) NegVector(v2,v2);     k<0 <=> from "inside" of object 
  166. */
  167. }
  168.  
  169.  
  170. /* This routine rotates a 2D point in a 2D system */
  171. /* I.e. (x,y) rotates (ar) radians counter clockwise around (0,0) */
  172.  
  173. void Rotate2D(double *x, double *y, double ar)
  174. {
  175.     double    x1,y1;
  176.  
  177.     x1=(*x)*cos(ar)-(*y)*sin(ar);
  178.     y1=(*x)*sin(ar)+(*y)*cos(ar);
  179.     *x=x1; *y=y1;
  180. }
  181.  
  182.  
  183. /* This routine rotates a 3D point in 3D space. The point (x,y,z) */
  184. /* is rotated around the x, the y and the z axis (in that order). */
  185. /* The vector RotV gives the degrees to rotate around each axis   */
  186.  
  187. void RotatePoint(POINT *p, VECTOR *RotV)
  188. {
  189.     double    x1,y1,z1;
  190.  
  191.     x1=p->x; y1=p->y; z1=p->z;
  192.     if(RotV->x!=0.0) Rotate2D(&y1,&z1,RotV->x);
  193.     if(RotV->y!=0.0) Rotate2D(&z1,&x1,RotV->y);
  194.     if(RotV->z!=0.0) Rotate2D(&x1,&y1,RotV->z);
  195.     p->x=x1; p->y=y1; p->z=z1;
  196. }
  197.  
  198.  
  199. /* This is the same thing, but with a vector instead of a point. */
  200.  
  201. void RotateVector(VECTOR *v, VECTOR *RotV)
  202. {
  203.     double    x1,y1,z1;
  204.  
  205.     x1=v->x; y1=v->y; z1=v->z;
  206.     if(RotV->x!=0.0) Rotate2D(&y1,&z1,RotV->x);
  207.     if(RotV->y!=0.0) Rotate2D(&z1,&x1,RotV->y);
  208.     if(RotV->z!=0.0) Rotate2D(&x1,&y1,RotV->z);
  209.     v->x=x1; v->y=y1; v->z=z1;
  210. }
  211.  
  212.  
  213. /* These routines will rotate in the reversed order (z,y,x). */
  214.  
  215. void RevRotatePoint(POINT *p, VECTOR *RotV)
  216. {
  217.     double    x1,y1,z1;
  218.  
  219.     x1=p->x; y1=p->y; z1=p->z;
  220.     if(RotV->z!=0.0) Rotate2D(&x1,&y1,RotV->z);
  221.     if(RotV->y!=0.0) Rotate2D(&z1,&x1,RotV->y);
  222.     if(RotV->x!=0.0) Rotate2D(&y1,&z1,RotV->x);
  223.     p->x=x1; p->y=y1; p->z=z1;
  224. }
  225.  
  226. void RevRotateVector(VECTOR *v, VECTOR *RotV)
  227. {
  228.     double    x1,y1,z1;
  229.  
  230.     x1=v->x; y1=v->y; z1=v->z;
  231.     if(RotV->z!=0.0) Rotate2D(&x1,&y1,RotV->z);
  232.     if(RotV->y!=0.0) Rotate2D(&z1,&x1,RotV->y);
  233.     if(RotV->x!=0.0) Rotate2D(&y1,&z1,RotV->x);
  234.     v->x=x1; v->y=y1; v->z=z1;
  235. }
  236.